import numpy as np
#https://blog.csdn.net/ReDreamme/article/details/108745255
#------------------平方根分解（Cholesky）分解#------------------
#计算下三角矩阵G-G=LD1/2，A=GGT-计算为方便示意不考虑时间空间代价
def Cholesky(matrix):
    w = matrix.shape[0]
    G = np.zeros((w,w))#实际上只用一半的空间就可以完成矩阵分解
    for i in range(w):
        G[i,i] = (matrix[i,i] - np.dot(G[i,:i],G[i,:i].T))**0.5
        for j in range(i+1,w):
            G[j,i] = (matrix[j,i] - np.dot(G[j,:i],G[i,:i].T))/G[i,i]
    return G

if __name__ == '__main__':
# A=np.random.randint(1,10,size=[3,3])  #注意A的顺序主子式大于零
    A = np.array([[4,2,1], [2,2,0], [1,0,3]] ) # 举一个例子
    G=Cholesky(A)
    print("原矩阵A：\n",A)
    print(G)
#------------------平方根分解（Cholesky）分解#------------------